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Abstract 

Mass-action kinetics is frequently used in systems biology to model 
the behaviour of interacting chemical species. Many important dynam- 
ical properties are known to hold for such systems if they are weakly 
reversible and have a low deficiency. In particular, the Deficiency Zero 
and Deficiency One Theorems guarantee strong regularity with regards 
to the number and stability of positive equilibrium states. It is also 
known that chemical reaction networks with disparate reaction struc- 
ture can exhibit the same qualitative dynamics. The theory of linear 
conjugacy encapsulates the cases where this relationship is captured by 
a linear transformation. In this paper, we propose a mixed-integer lin- 
ear programming algorithm capable of determining weakly reversible 
reaction networks with a minimal deficiency which are linearly conju- 
gate to a given reaction network. 

Keywords: chemical kinetics; stability theory; weak reversibility; linear 

programming; dynamical equivalence 

AMS Subject Classifications: 80A30, 90C35. 



1 Introduction 

A chemical reaction network is given by sets of chemical reactants reacting 
at prescribed rates to form sets of chemical products. Under suitable as- 
sumptions, such as mass-action kinetics and spatial homogeneity, the time 
evolution of the concentrations of the chemical species can be modeled by 
a set of autonomous polynomial ordinary differential equations. Such mass- 
action systems are frequently used to model systems in systems biology and 



other areas of computational biology 27,35, 36 



The systematic study of chemical reaction networks was initiated in 1972 
in the papers [7, 16, 18 . In 18 , the authors presented a condition on the 



positive equilibrium concentrations, called complex balancing, which is suf- 
ficient to guarantee that within each linear invariant space of the network 
there exists a unique equilibrium concentration, and that this concentration 
is locally asymptotically stable relative to its invariant space. In [T] and [l6] , 
the authors related the capacity of a network to exhibit complex balancing to 
a nonnegative network parameter called the deficiency. In particular, they 
showed that a network is complex balanced for all sets of rate constant values 



if and only if it is weakly reversible and has a deficiency of zero. This defi- 
ciency-onented approach to analysing chemical reaction networks has since 
been applied to a wide variety of biochemical systems, including enzymatic 
models, signal transduction, and phosphorylation networks [4| |24[[25|[30l i [3T] . 
It is also known that many qualitative properties of the dynamics of reac- 
tion networks can be shared by networks with disparate network structure. 
A thorough study of dynamical equivalence, the property that two reaction 
networks give rise to identical mass-action dynamics, was conducted in L2j. 
In [21], dynamical equivalence was extended to linear conjugacy, whereby 
the trajectories of two mass-action systems could be related by a non-trivial 
linear transformation. The problem of determining dynamically equivalent 
networks with the greatest and fewest number of reactions was placed in a 



mixed-integer linear programming (MILP) context in 32 . This methodol- 
ogy has been extended to linear conjugacy in the papers [23] and 22 . 

Linear conjugacy effectively creates classes of mechanisms which, de- 
spite disparate network structure and properties, have equivalent dynamics. 
A long-standing problem in chemical reaction network theory, first stated 
in 15 , is to "... look for a mechanism in a class of mechanisms with a 
given - chemically relevant - property. Such a property may be conservativ- 
ity, (weak) reversibility, zero deficiency or just structural stability as well." 
In this paper, we use linear conjugacy and MILP methods to answer this 
challenge for the class of weakly reversible mechanisms where the structural 
property of interest is the deficiency. 

In general, we will be interested in mechanisms within the class of lin- 
early conjugate networks with minimal deficiency. This is due to the ob- 
servation that the behaviour of a chemical reaction network's mass-action 
system is generally more regular (e.g. fewer steady states, lower capacity 
for oscillations, etc.) for networks with a lower deficiency than a higher 
deficiency iTl^flTlflG] . Consequently, within the class of mechanisms which 
are linearly conjugate, if there are mechanisms with differing deficiencies, 
the network with the lowest deficiency is generally preferred. We therefore 
modify the existing MILP framework to compute weakly reversible linearly 
conjugate networks which have a minimal deficiency. The methodology is 
based on known properties of the kernel of the kinetics matrix of a weakly 
reversible chemical reaction network. 



2 Background 

In this section we present terminology and notation relevant to the study of 
chemical reaction networks. We introduce the notion of the deficiency of a 
network and a few classical results which relate the deficiency of a network 
to the dynamics of the network's corresponding mass-action system. We 
also introduce the notion of two chemical reaction networks being linearly 
conjugate. 

2.1 Chemical Reaction Networks 

The chemical species or reactants of a network will be given by the set 
S = {Xi, X2, . . . , Xn}- The combined elements on the left-hand and right- 
hand side of a reaction are given by linear combinations of these species. 
These combined terms are called complexes and will be denoted by the set 
C = {Ci, C2, . . . , Cm} where 

n 

Ci = '^aijXj, i = l,...,m 
i=i 

and the aij are nonnegative integers called the stoichiometric coefficients. 
The complex with all stoichiometric coefficients equal to zero will be called 
the null complex and denoted by Ci = 0. We define the reaction set to 
be 7?. = {{Ci,Cj) I Ci reacts to form Cj}. The property {Ci,Cj) G TZ will 
commonly be denoted Cj — ;■ Cj. To each {Ci,Cj) G 7^ we will associate 
a positive rate constant k{i,j) > and to each {Ci,Cj) ^ TZ we will set 
k{i,j) = 0. The triplet M = {S,C,TZ) will be called the chemical reaction 
network. 

The above formulation naturally gives rise to a directed graph G{V, E) 
where the set of vertices is given hy V = C, the set of directed edges is given 
hy E = TZ, and the rate constants k(i,j) correspond to the weights of the 
edges from d to Cj. In the literature this has been termed the reaction 
graph of the network |18] . Since complexes may be involved in more than 
one reaction, as a product or a reactant, there is further graph theory we 
may consider. A linkage class is a maximally connected set of complexes, 
that is to say, two complexes are in the same linkage class if and only if 
there is a sequence of reactions in the reaction graph (of either direction) 
which connects them. We will denote by i the number of linkage classes 
in a network. A reaction network is called weakly reversible if Cj — )• Cj 
for any Ci,Cj £ C implies there is some sequence of complexes such that 



A directed graph is called strongly connected if there exists a directed 
path from each vertex to every other vertex. A strongly connected component 
of a directed graph is a maximal set of vertices for which paths exist from 
each vertex in the set to every other vertex in the set. A strongly connected 
component is called terminal if there is no reaction leading from a vertex in 
the strongly connected component to a vertex not in the component. For a 
weakly reversible network, all strongly connected components are terminal, 
and they correspond to the linkage classes of the reaction graph. 

Assuming mass-action kinetics, the dynamics of the specie concentra- 
tions over time is governed by the set of differential equations 

f = >--A.-*(x) (1) 

where x = [xi 2:2 ■ ■ ■ ^n]'^ is the vector of reactant concentrations. The 
stoichiometric matrix Y contains entries [Y]ij = Uji and the Kirchhoff or 
kinetics matrix Ai^ is given by 

I A;(J,^) if z^j 

for i,j = l,...,m. When we speak of the structure of a kinetics matrix, we 
are referring to the distribution of positive and zero entries, which determines 
the network structure of the corresponding reaction graph. Finally, the 
mass-action vector ^(x) is given by 

n 

^.•(x)=n^p' j=i'---'"^- (3) 

i=l 

The behaviour of solutions of M can be further understood by consid- 
eration of the reaction vectors 

[y].,j- - [y].,, for iCi,Cj) £71 

otherwise 

where [y].,i denotes the ith column of Y. The span of the vectors is called 
the stoichiometric subspace and is denoted by 5* = span{vij \ {Ci, Cj) E TZ}. 
We will denote the dimension of 5 by s = dim(S'). 

It is clear that the right-hand side of M is contained in S so that the vec- 
tor field always directs trajectories within an affine translation of S. It can 
be shown that trajectories are restricted to the stoichiometric compatibility 



classes (xq + S) n R^q 18 , 37 



2.2 Complex Balanced Networks and Deficiency 

The structure of the reaction graph of a chemical reaction network often 
plays an important role in determining the dynamical behaviour of trajec- 
tories of ([I]). 

A particularly important structural parameter of a chemical reaction 
network is the deficiency, which was introduced in [7, 16 and has been 
studied extensively since (§[9||l2||29] . 

Definition 2.1. The deficiency of a chemical reaction network is given by 

5 = m — i — s 

where m is the number of stoichiometrically distinct complexes, i is the num- 
ber of linkage classes, and s is the dimension of the stoichiometric subspace. 

It is easy to show that the deficiency of a network may only take on non- 
negative integer values. 

The deficiency of a network is often related to properties of the equilib- 
rium set of the mass-action system (II| . It is strongly related to the capacity 
for a network to admit complex balanced equilibrium concentrations. 

Definition 2.2. An equilibrium concentration x* G R"q is called a complex 
balanced equilibrium concentration if 

Ak • ^(x*) = 0. 

A chemical reaction network is called complex balanced if every equilib- 
rium concentration is a complex balanced equilibrium concentration. 

It is known that if a network is complex balanced at one equilibrium con- 
centration then it is complex balanced at all of them (Theorem 6A, il8i)- 
Consequently, a network for which any equilibrium concentration is complex 
balanced is a complex balanced network. 

The following results relate the deficiency of a network to its capacity 
for complex balancing (see Theorem 4A of [16| and Theorem 4.1 of |7j). 



Theorem 2.1 (Deficiency Zero Theorem). A chemical reaction network 
is complex balanced for all rate constant values if and only if it is weakly 
reversible and has a deficiency of zero. 

Theorem 2.2. // a mass-action system is weakly reversible, then the defi- 
ciency corresponds to the number of algebraically independent conditions on 
the rate constants which need to be satisfied in order for the system to be 
complex balanced. 



Theorem 2.1 is a special case of Theorem 2.2, taking the deficiency to be 
zero. We state is separately, however, due to its historical importance. 

Strong dynamical properties are known to follow from complex balanc- 
ing; in particular, the following was proved in jl8|. 



Theorem 2.3 (Theorem 6A and Lemma 4C, [18] )• If a chemical reaction 
network is complex balanced then there exists within each positive stoichio- 
metric compatibility class exactly one equilibrium concentration, and that 
equilibrium concentration is locally asymptotically stable relative to its com- 
patibility class. 



The surprising implication of Theorem 2.1 and Theorem 2.3 is that we 
can know very strong properties about the equilibrium set of the mass-action 
system ([I]) based on structural properties of the reaction graph alone. Most 
strikingly, the results hold regardless of the values of the rate constants. 
This is particularly important for biological examples, where the relevant 
rate constants are frequently unknown or unknowable. In cases where the 
deficiency in nonzero, supplemental conditions on the rate constants are 
required to bridge the gap from weak reversibility to complex balancing. 
Theorem |2.2| tells us that we are more likely to be interested in networks 
with a lower deficiency, since networks with a lower deficiency are, in some 
senses, closer to complex balancing (and the regular dynamics guaranteed 



by Theorem 2.3) than networks with a higher deficiency. 

A result which further relates deficiency to properties of the equilibrium 
set of a chemical reaction network is the following, which can be found 
in 11 . It is worth noting, once again, that a lower deficiency generally 
guarantees more predictable behaviour than a higher deficiency. (Further 
work investigating the capacity of a deficiency one chemical reaction network 
to exhibit multistability was presented in jol.) 

Theorem 2.4 (Deficiency One Theorem). Consider a chemical reaction 
network with deficiency 5 and linkage classes Ci, . . . , Ci. Let 6i, i = 1, . . . ,i, 
be the deficiency of the linkage class Ci considered as its own subnetwork. 
Suppose that 

1- Si <1, fori = l,...,£; 

^- ^i=i Si = S; and 

3. each linkage class contains only one terminal strongly connected com- 
ponent. 



Then, if the network admits a positive equilibrium concentration, there is 
exactly one equilibrium concentration in each positive stoichiometric com- 
patibility class. Furthermore, if the network is weakly reversible, then for 
every set of rate constants the network admits a positive equilibrium concen- 
tration. 

2.3 Linearly Conjugate Networks 

Under the assumption of mass-action kinetics, it is possible for the trajecto- 
ries of two reaction networks M and A/"' to be related by a linear transforma- 
tion and therefore share many of the same qualitative features (e.g. number 
and stability of equilibria) . This phenomenon was termed linear conjugacy 
in [21]. 

For completeness, we include the formal definition of linear conjugacy as 
presented in ^21j. We will let $(xo,t) denote the flow of ([I]) associated with 
7V^ and ^(xq, t) denote the flow of ([l]) associated with A/"'. 

Definition 2.3. We will say two chemical reaction networks N andN' are 
linearly conjugate if there exists an infective and surjective linear mapping 
h : W^Q ^ RIq such that h($(xo,t)) = ^(h(xo),t) for all xq G M^^q- 

It is known that an injective and surjective linear transformations h : ]R"q i— )• 
M"g can consist of at most positive scaling and reindexing of coordinates 
(Lemma 3.1, l21j). Linear conjugacy has been subsequently studied from a 
computational point of view in |23| and |22|. 

Linear conjugacy is a generalization of the concept of dynamical equiv- 
alence whereby two reaction networks with different topological network 
structure can generate the same exact set of differential equations (IT]) . Two 
dynamically equivalent networks M and M' are said to be alternative real- 
izations of the kinetics M , although it is sometimes preferable to say that 
M is an alternative realization of M' or vice- versa. Since the case of two 
networks being realizations of the same kinetics is encompassed as a special 
case of linear conjugacy taking the transformation to be the identity, we will 
focus on linearly conjugate networks. 

If a network M can be shown to be linearly conjugate to a network J\f' 
from a class of networks with desired properties, those properties of M' 
are transferred to M. This raises the question of how to find a linearly 
conjugate network A/"' when only the original network M is given. This was 



studied in 23 and |22l where the authors built upon the linear programming 



algorithm introduced in [32]. We can impose that a network J\f' be linearly 



conjugate to a given network N with the set of Unear constraints 



(LC) { 



^[A]ij=0, j = l,...,m 

< [Ablij < 1/e, i,j = l,...,m, i^ j 
[Af,]ii < 0, i = 1,. . .,m 
1, 



(4) 



e<Cj<l/e, j = i,...,n 

where < e <C 1, and the matrices M e M^x"* and T G M"^'^ are given by: 

M = y-^fc, and (5) 

T = diagjc}. (6) 

The kinetics matrix for the network J\f' can by constructed from Ai, G ]g^"^^^ 
and c G R"q by the relation 

A', = Ab-dmg{^{c)}. (7) 

Finding a network satisfying Ml) and then solving ([7| is sufficient to 
determine a network A/"' which is linearly conjugate to M via the transfor- 
mation h(x) = T~^x. It is important to note that, while the matrix A^ does 
not exactly correspond to the kinetics matrix A^ of the network M' , it is 
a kinetics matrix and has the same structure as A'f,. Most importantly, Ai, 
corresponds to a weakly reversible network if and only if the A^ defined by 
^ does. 

It is also worth noting that, without further clarification, the value m 
may be different in the linear conjugacy constraints Q than that in Def- 
inition 2.1 The value m in the constraints Q corresponds to the set of 
potential complexes represented in the matrix Y. This set must be fixed 
before the optimization procedure can begin and as such there is no way to 
guarantee a priori that all of the complexes assigned to Y will be present in 
the network M' . The value may consequently be different than the value m 



used in Definition 2.1 to represent the number of complexes in the network 
AA'. (It is common to initialize Y with the reactant and product complexes 



from the original network A^ [23 32-34 ; however this choice is by no means 
necessarily the best [22|.) 

We will avoid this notational difficulty by expanding the definition of 
the network J\f' to include all of the complexes initialized in Y; that is to 
say, we will allow complexes to be in A/"' even if they are not the reactant or 
product complex for any reaction in J\f' . Consequently, the values of m in 
the linear constraints and the value of m in Definition 12.11 will coincide for 
the network M' . 



3 Minimal Deficiency Networks 

We saw in Section [2.2| that many properties of the equilibrium set of a mass- 
action system (II| can be related to the value of the deficiency of the cor- 



responding chemical reaction network. Theorem |2 . 1| and Theorem 2.2 show 



that weakly reversible networks with lower deficiency are closer to complex 



balanced networks (and the regular dynamics guaranteed by Theorem 2.3) 



than ones with higher deficiency. Theorem 2.4 gives conditions under which 
a network with a nonzero deficiency can still be guaranteed to have a unique 
positive equilibrium concentration in each compatibility class. 

In both of these cases, we find that for weakly reversible networks a lower 
deficiency is more likely to be related to regular dynamics than a higher de- 
ficiency. In cases where a network is linearly conjugate to multiple weakly 
reversible networks, therefore, we are likely to be most interested in the net- 
work with the minimal deficiency. In this section, we develop a mixed-integer 
linear programming procedure capable of determining a linearly conjugate 
weakly reversible network with a minimal deficiency. 

3.1 Parameters of Interest 



We recall that, according to Definition 2.1, the deficiency of a chemical 
reaction network depends on three structural parameters of the reaction 
network: the number of stoichiometrically distinct complexes m, the number 
of linkage classes i, and the dimension of the stoichiometric space s. 

In order to avoid ambiguity with the value m used to represent the 
number of potential complexes in the linear constraints Q and the value in 



Definition 2.1 we allow A/"' to contain complexes which do not correspond 
to the reactant or product complex of a reaction in A/"'. These unused 
complexes will appear in the reaction graph as isolated nodes. As such, 
each unused complex will correspond to a linkage class in itself and this 
linkage class will be trivially strongly connected. Consequently, including 
these unused complexes in A/"' will not change the value of the deficiency of 
M' , since the increase to m will be offset by a corresponding increase in i 



in Definition 2.1 Including these complexes will also not alter properties 
related to the reversibility of the network. 

In order to minimize the deficiency within a linear programming frame- 
work, we need to find a way to quantify the parameters m, i, and s. We 
notice that, since linear conjugacy preserves the dimension of invariant lin- 
ear spaces, the dimension of the stoichiometric space s is fixed for M' by 
the dimension of the largest invariant linear space of the original network 



10 



N . Since N' is weakly reversible, the parameter s corresponds to the di- 
mension of the stoichiometric space (see Corollary 1 of |17|); however, this 
correspondence is not necessarily true of AA (see Example 2 of [19|). 

We also notice that the value of m is predetermined for the optimization 
procedure. Consequently, in order to minimize the deficiency it is sufficient 
to maximize L In other words, we need to maximize the number of linkage 
classes, where each unused complex corresponds trivially to its own linkage 
class. 

3.2 Counting Linkage Classes 

We need to keep track of the number of linkage classes. In general, this is a 
difficult task; however, we are aided by the following result. 



Theorem 3.1 (Theorem 3.1 of [14] and Proposition 4.1 of |6|). Let A^ be a 

kinetics matrix and let Ai, i = 1, . . . ,£, denote the support of the i*'' linkage 
class. Then the reaction graph corresponding to Ak is weakly reversible if and 
only if there is a basis ofker{Ak), (b^^), . . . ,b(^-'}, such that, fori = !,...,£, 

\ bf =0, jf A,. 

It is easy to see that, for a weakly reversible network, the dimension of 



ker(^fc) given by Theorem 3.1 corresponds to the number of linkage classes. 



Consequently, the parameter £ here coincides with the earlier usage. 



It is typical to assume in Theorem 3.1 that the kinetics matrix A^ cor- 
responds to a network where every complex appears on either the reactant 
or product side of at least one reaction. It is easy to extend this to the 
case where complexes are not used by the reaction network by noting that 
any unused complexes will contribute an element to ker(^fc) satisfying ([s]) 
corresponding to a single positive value in the coordinate corresponding to 
the unused complex and zeroes elsewhere. In other words, for weakly re- 
versible networks with unused complexes, we can extend the basis of ker(^fc) 
by considering unused complexes as their own linkage classes. Throughout 
this section, we will allow i to correspond to both traditional linkage classes 
containing several complexes and unused complexes. 



Theorem 3.1 implies that for a weakly reversible network, the supports 



of the basis elements of ker(^fc) represent a complete partition of the set 



11 



{1, . . . , m}, that is to say, we require 

Afci n Afca = 0, for all fci, /C2 = 1, • • • , i, h / /c2 

[jAk = {l,...,m}. ^^) 

fc=i 

The value of i, however, is not known; in fact, it is what is to be de- 
termined through the procedure. In order to implement this into a compu- 
tational framework, therefore, we need to determine an upper limit for the 
number of possible supports A^ . We recall that the deficiency is a nonnega- 
tive parameter, so that we have 6 = m — i — s>0 [7|[l6]. Consequently, we 
have i < m — s and therefore may use m, — s as an upper bound for i. It is 
also clear that the deficiency will be zero if and only if this upper bound is 
attained. 

We now introduce the binary variables 7jfc G {0,1}, for i = l,...,m, 
k = 1, . . . ,m — s, defined according to 

/ 1, if i G Afe , . 

^^'= = 10, if.^A,. (10) 

The variables ■jik keep track of how the supports of the basis vectors in 
ker(^fc) according to ([s]) partition the set {1, . . . , m}. For weakly reversible 
networks, this corresponds to an assignment between the complexes and 



the linkage classes by Theorem 3.1 In other words, jik = 1 if and only if 
Ci € Ck- 

We also introduce variables 6k € [0,1], k = 1, ... , m—s, defined according 

to 

1, ifsupp(Afc)/0 

0, if supp(Afc) =0. ^ ^ 

The variables 6k keep track of whether the A:th partition of { 1 , . . . , m} is 
empty or nonempty. It should be noted that, while we would like the 6k s 
to count the number of non-empty supports, and are therefore interested in 
only the values 6k = and 6k = I, it will be possible to relax the integrality 
of the 6k s to vary continuously within the range [0, 1]. This will be justified 
in Section [3?5l 

In order to accommodate the complete partition requirements ([9| as 



12 



linear constraints, and to accommodate (11), we impose 



(CP) { 



,m 



m 



(12) 



m 



^ 7ifc = 1, i = l,..., 
fc=i 

m 

^7ifc-e6'fc>0, k = l, 

i=l 

m ^ 

-^lik + -Ok>Q, fc = l,..., 

k=l 
7jfcG{0, 1}, i = 1,... ,m,k = 1,. . .,m 

6'fce[0, 1], k = l,...,m-s. 



where < e ^ 1 is sufficiently small and can be chosen to be the same as 
The first constraint set guarantees that each complex appears in exactly 



one partition. The second two constraint sets of (12) correspond to the 
constraint 



IIL ^ 

o<eek<Y,^ik<-^ 



i=l 



which keeps track of whether the /cth partition is empty or nonempty. If 
no element is in the kth partition, then the sum is zero, which forces 9^ to 
be zero (first inequality). If there is an element in the fcth partition, then 
the sum is nonzero, which forces 6^ to be nonzero (second inequality). An 



argument in Section 3.5 will allow us to conclude that any nonzero 9k must 



be one, so that this fulfills the requirements of (11). 



3.3 Constructing the Kernel 



We still need to guarantee that the sets A^ considered in Section 3.2 



cor- 



respond to the supports of vectors in ker(j4fe) which satisfy ([8]). In other 
words, we need to restrict ourselves to sets Aj, i = l,...,i, where there 
exists a b^*) satisfing m« and 



A, ■ b« = 0, 



1, 



(13) 



We follow the technique outlined in the paper |23| for determining weakly 
reversible networks. We define a matrix $ G j^mxm, ^j^j^ entries 



^.. 



[Mij ■ bj 



(14) 



13 



where b = X]^,=x b*^" and {b'^-*, . . . , b^^} is any set of vectors satisfying ish 
and forming a basis of ker{Af,). We can see that the system of non-hnear 
equations p^ is satisfied if and only if the system of hnear equations 



y ^ ^ij = 0, for ah z = 1, . . . , m, A; = 1 



, . . . , -u, 



(15) 



is satisfied. 

We know that $ is a kinetics matrix since A), is and (14) preserves this 
property. That is to say, we have 



V^ $ji = 0, for ah i = 1, . . . , 



m 



(16) 



and 



$jj > 0, for alH, j = 1, . . . , m, z 7^ j. 



Consequently, we can solve the diagonal elements of $ in ( 16 ) and substitute 



them into (15) to get the simplified constraint set 



y^ ^ij = ^ ^ji, for i = 1, . . . , m, A; = 1, . . . 



(17) 



3+i 






We need to derive linear constraints capable of constructing a matrix $ 



according to (14) which satisfies ( 17) and has the same structure as Ah. This 



is made more challenging than the case considered in |23j by the requirement 
that the kernel vector b in (14) decompose according to the partitions A^, 
k = !,...,£. 

We notice first of all that we do not know how many partitions A^ we will 
need, so we will take the upper bound m — s on the number of partitions (see 
Section 3.2). We also notice that the construction (14) requires that $jj = 
for any i and j such that i S A^ and j A^ for some k = 1, . . . ,m — s. 
That is to say, we do not permit reactions to occur between complexes in 
different linkage classes. This can be accommodated by noticing that, if two 
indices i and j are on the same support Afc, then we will have 7^^ — ^jk = 
for all k = 1, . . . ,m — s, where as if i and j are not on the same support 
Afc, we will have ^ik — 'Jjk G {~1) 0, 1} and attain the value —1 for at least 
one k. Consequently, we can accommodate this requirement with the linear 
constraints 



'^ij < l/e(7ifc -7jfc + 1), 
ioT i,j = l,...,m,iy^ j,k = l,... 



,m 



(18) 



14 



We can accommodate (17), (18), and the requirement that $ have the 



same structure as Ah, with the Hnear constraint set 



(Ker) < 



E^^' = E^' 






1 



1=1 



^ij < -{lik -ijk + 1) 

^ij > e[Ab]ij 
I, 



(19) 



^^J< 



\A 



b\ij 



hJ 



1, 



<m,iy^j, k = l, 



,m — s. 



Notice that (17) can be generahzed to the first constraint in (19) because of 



the imposition that ^u = for any Z ^ A^ in the left sum. This is guaranteed 



by (18). 



3.4 Uniqueness of Solution 

The constraint sets (B, (12), and (19), form the basis of a mixed-integer 
linear programming (MILP) problem. This is due to the non-continuous 
binary variables ^ik, i = 1, . . . ,m, j = 1, . . . ,m — s, which keep track of how 
the complexes are assigned to the partitions. 

MILP problems are known to NP-hard and are generally approached by 
a branch-and-bound method (for an accessible introduction to branch-and- 
bound methodology, see |13[|26| ). One well-known complicating factor for 
branch-and-bound methods is non-uniqueness of the integer portion of the 
problem. There is nothing in the constraint sets (12) which guarantee a 
unique assignment of the complexes to partitions. That is to say, a link- 
age class could be assigned to the first partition as easily as the second or 
third. Consequently, we would like to introduce further constraints which 
guarantee a unique partitioning structure for the optimal solution. 

An intuitive way to structure the partitioning variables jik is to always 
assign the first complex to the first partition, and then to assign each sub- 
sequent free complex to the next available partition. That is to say, if the 
second complex is not in the same partition as the first complex, it will be 
assigned to the second partition. If it is in the same partition as the first 
complexes, but the third is not, then the third complex will be assigned to 
the second partition. 
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This unique partition structure can be guaranteed by the Hnear con- 
straint set 



(U) 



i— 1 m~s 

j=l l=k+l 

i = 1, . . . ,m,k = 1, . . . ,m — s,k < i. 



(20) 



This constraint set guarantees that, in ascending order, each complex 
which is assigned to a new partition is assigned to the unused partition with 
the lowest index. It does so in the following way: 



1. 



For the case i = k = 1, the left sum in (20) is empty and therefore 
returns a value of zero. This forces 71; = for I = 2, . . . ,m — s, which 
forces 711 = 1 by the first condition of (12). In other words, the first 



complex is always assigned to the first partition. 



2. The left sum in (20) is zero for i = k = 2, which forces 72^ = for 
3, . . . ,m — s. This forces the second complex to be assigned to 



J 

either the first partition or the second partition by (12), depending on 



whether it is grouped with the first complex or not. 
3. Following an induction on the complexes as i = 3, . . . ,m, we can see 



that the constraint (20) corresponding to the ith. complex and the 



first unused partition will always guarantee this complex may not be 
assigned to a partition of a higher index than the first unused partition. 
Consequently, if it not grouped with an earlier complex, it must be 
assigned the index of the next available partition. 



4. It is easy to see that the conditions (20) are satisfied for all of the 



entries not corresponding to leading ones, so that we are done. 



3.5 Minimizing the Deficiency 

We know that minimizing the deficiency 6 is equivalent to maximizing the 
number of linkage classes i, where we include the unused complexes as link- 
age classes unto themselves. 

3.2| we introduced variables Oi^. € [0, 1] , A; = 1, . . . , m — s, ac- 

,m} 



In Section 



cording to (ITl ) to keep track of the number of partitions A^ of { 1 , 
which corresponded to vectors in ker(^i,) satisfying m^. We have no guar- 
antee, however, that these sets correspond to a complete basis of ker(^;,) 
since any two vectors in ker(^fe) satisfying ([s]) can be added to one another 
to produce another vector in ker{Ab) also satisfying ([8]). (This corresponds 
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to multiple linkage classes being placed on the same support A^.) We also 
have no guarantee that the ^^'s will properly enumerate the number of basis 
elements of ker(Ab) since any 6^ corresponding to a nonempty support may 
attain a value anywhere between zero and one, and not just the value of 
one. 

We need to guarantee that the sum of the ^^'s corresponds to the max- 
imal partition of {1, . . . ,m} and that each 9k attains a value of one when 
it corresponds to a nonempty support. Consider the objective function cor- 
responding to maximizing the number of linkage classes, which corresponds 
to minimizing the deficiency. We have 

{m— 1 
minimize ^ -Ok (21) 

fc=i 

over the constraint sets Q, (12), and (19). We notice that this objective 
function guarantees that, if a 0^ is in between zero and one, it will attain 
the value one and that a maximal partition of {1, . . . ,m} will be chosen, 
because in both cases such a situation is more optimal than the alternative. 
In other words, the Ok's properly count the number of linkage classes when 



(21) is imposed. 

If the off-diagonal elements are solved for in ([4]) , the algorithm presented 
here for finding a weakly reversible network with minimal deficiency which is 
linearly conjugate to a given network contains m{2m — l)+n — s continuous 
decision variables and m{m — s) binary decisions variables. 

It should be noted that a linearly conjugate network with the maximal 
deficiency cannot be obtained in the same fashion as outlined here. This is 



due to having to maximize the sum in (21 ) in order to make the correspon- 
dence between the sum of the 0fc's and the value of £. Without this step, we 
could not rule out linkage classes clumping together onto the same support 
in ker(^fc). 

4 Examples 

We now introduce a few examples which illustrate the methods presented 
in this paper. All computations are performed on the primary author's 
personal-use Acer laptop (AMD Athlon II Neo K125 Processor 1.70 GHz, 4 
GB RAM). 

Example 1: We propose the following mechanism. Consider a substrate 
which can bind to an enzyme T at one of three binding sites and let Tioo, 
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Tolo, and Tqoi denote the enzyme with the substrate bound at the first, 
second, and third binding site, respectively. Suppose that binary cohisions 
between the substrates can cause a spontaneous shift in the substrate from 
one binding site to another, but no transfer of substrate from one bound 
enzyme to another. This web of interactions can be visuahzed by the network 
in Figure [l] 



2T 

^ ' 100 

H % 

Tioo+Tqo X ^ '100+ '010 

H % H % 



' OIO"*" ' 001 ^ 010 



Figure 1: Network of singularly bound enzymes with three binding sites. In- 
teractions between the enzymes allow transfer of substrates from one binding 
site to another. 

This network is weakly reversible and has a deficiency of three {jn = 
6,£ = 1, s = 2, (5 = m — £ — s = 3). Due to the high deficiency, we 



may not apply Theorem 2.1 or Theorem |2.4[ however, since the network is 



weakly reversible, we may apply Theorem 2.2 This tells us that there are 
three algebraically independent conditions on the rate constants which must 
be satisfied in order for the system to be complex balanced and therefore 
fall within the scope of Theorem 2^ If we set Ci = 2Tioo, C2 = 2Toio, 
C3 = 2rooi, C4 = Tioo + 7bio> C's = ^100 + ^001) and Cg = Tqio + Tboi, the 
required conditions are K1K2 = K^, K1K3 = K^, and K2K3 = Kq, where 
Ki is the absolute value of the {i x i) minor of A^. (For details, see 111.) 

We might wonder what additional behaviour is permitted by the reaction 
network given in Figure [T} A general analysis can be conducted by the CRN 
toolbox made freely available online 20 . This is a powerful computational 



toolbox capable of determining whether a chemical reaction network has 
the capacity for zero eigenvalues, injectivity, multistability, and concordance 
(5|[8| |10||28| . The toolbox's Higher Deficiency Report reveals that the mass- 
action system corresponding to the network has the capacity for multiple 
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positive steady states for the rate constant choices 



k{l,4) = 7.389056 k{4, 1) = 2.7182818 k{5, 4) = 4.3002585 

k{l, 5) = 7.389056 k{4, 2) = 2.7182818 A;(5, 6) = 4.3002585 

k{2, 4) = 1 /c(4, 5) = 55.125832 A;(6, 2) = 32.08195 

k{2, 6) = 1 /c(4, 6) = 29.019118 A:(6, 3) = 1.5819767 

/c(3, 5) = 2.5026503 A;(5, 1) = 45.90757 A;(6, 4) = 1.5819767 

k{3, 6) = 2.5026503 A;(5, 3) = 4.3002585 A;(6, 5) = 1.5819767. 



(22) 



For rate constant choices lying away from these sets, the dynamics is 
significantly more difficult to analyze. We consider the deficiency-reducing 
methodology introduced in this paper on two sets of rate constant choices 



(in addition to (22)). The two additional sets of rate constant choices are 



kii,j)=j, «,J = l,---6, ii,j)eTZ 



and 



Hhj) 



^,J = l,---,6, {i,j)£n. 



(23) 



(24) 



In other words, we consider the rate constant choices where all the rates 
going into a particular complex are the same, and the rate constant choices 
where all the rates going out of a particular complex are the same. For 
simplicity, we take the differences between complexes to scale according to 
the index of the complex itself. 

We optimize (21) over the contraint sets (Q, ( [T2| ), (19), and (20), in 
GLPK. The algorithm runs in under a second for all rate constant sets. 



The rate constants choice (22) produces a network with deficiency three, 



and consequently we conclude that there is no network which is linearly 
conjugate to the network in Figure [T] with a lower deficiency. 

The results of the optimization for ( |23[ ) and ( 24 ) are contained in Figure 
[2] In both cases, the conjugacy constants are ci = C2 = C3 = 1. The 



network corresponding to the rate constant set ( 23 ) is deficiency two while 



the network corresponding to the rate constant set (24) is deficiency one. 



Neither network is amenable to application of Theorem |2.1| or Theorem 



2.4, The CRN toolbox can, however, by used to verify (by the algorithm 
presented in [9j ) that the network contained in Figure W[h) admits at most 
one equilibrium concentration in each positive compatibility class. (Since 
weakly reversible networks contain at least one equilibrium concentration in 
each compatibility class, this is sufficient to guarantee that the network has 
exactly one equilibrium concentration in each compatibility class [3].) 
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2T 



100 



<^> n \ 



2T 

^ ' 001 " 
' 100+ ' 001' 



2T 



010 



' 100+ ' 010 



' 010+ ' I 



^ ' ^ ' 100 <C T010+T001 



^ ' 010 <[ ^ ' 100+ ' 001 

^ ' 001^ ' 100+ ' 010 



001 



Figure 2: Weakly reversible networks which are dynamically equivalent to 
the one contained in Figure fl} The network in (a) is dynamically equivalent 
for the rate constant choices (|23|) and has a deficiency of two. The network 



in (b) is dynamically equivalent for the rate constant choices (24) and has 
a deficiency of one. 



The rate constant values corresponding to the two networks in Figure [2] 



are 



(a) A:(l,2) 
fe(l,3) 
fc(2,l) 
fe(2,3) 



2 A;(3,l) 

2.5 A;(3,2) 

2 A;(4,5) 

3 fe(4,6) 



2.5 /c(5,4) = 2 

3 /c(5,6) = 3 

4 A;(6,4) = 3 
7 A:(6,5) = 6, 



(5) /i:(l,6) = l A;(3,4) = 3 A;(5,2) = 5 
/c(2,5) = 2 A;(4,3) = 4 A;(6, 1) = 6. 

Example 2: Consider the set of differential equations given by 



dX2 

~dF 

dX3 

dt 



1- 

2x1 



Xi 



Xl + X2X3 



2x2X3 - 2x2 + 2x1 



Xl - X2X3 + X2 



X3. 



(25) 



We can produce a chemical reaction network which generates ( 25 ) using the 



algorithm presented in ll5i and reproduced in [33i. This gives a network 



20 



involving the complexes 



Ci = 0, (^2 = Xi, C3 = 2X1, C4 = X2 + X3, 
C^ = Xi + X2 + X3, Ce = Xi + X2, C7 = X^, 
Cg = 2X2, Cg = X2, Cio = 2X3, Cii = X2 + 2X3, 
C12 = ^1 + -^^35 Cl3 = 2X2 + X3. 

The corresponding network is not weakly reversible and has a deficiency 
of eight. It is therefore not amenable to Theorem 2.2 or Theorem 2.4 



In order to apply these results, therefore, we would like to find a weakly 
reversible network which is linearly conjugate to this one and has a minimal 
deficiency. We construct the matrices 



Y 



0120110000010 
0001110210102 
0001101002211 



and 



M 



1 


-1 


-1 


10 

















2 





-2000 


-2 





2 








1 





-10 


1 





-1 






and optimize (21) over the contraint sets Q, (12), (19), and (20). The 
upper limit for the number of partitions, m — s, can be easily determined to 
be 10. 

A quick computation in GLPK produces the network given in Figure 
p[a). We can readily see that this is a weakly reversible zero deficiency 
system and therefore falls within the scope of the networks considered by 



Theorem 2.1 and Theorem 2.3 Consequently, we know that (25) has exactly 



one positive equilibrium concentration and that this equilibrium concentra- 
tion is locally asymptotically stable. 

It is interesting to note that the network contained in Figure p[a) is not 
the only weakly reversible network capable of producing (25). If we do not 



insist on maximizing the number of linkage classes, other networks can be 
selected with a sub-optimal deficiency value. The network in Figure [3^ 



also linearly conjugate to a network generating (25 ) with the same conjugacy 



IS 



constants as that of Figure pf a). Despite being weakly reversible, however, 
this network has a deficiency of two and is not amenable to either Theorem 
Oor Theorem [231 
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x,^x,+X3 x,:^x,+x,^ 

2X J=i2X '■^ '^ 

"^^2 <-;^ ^^' 2X, 2X, 



Figure 3: Two weakly reversible network structures which are linearly con- 



jugate to a network generating kinetics (25). The conjugacy constants are 
ci = C3 = 1 and C2 = 2. The network in part (a) has a deficiency of zero 
while the network in part (b) has a deficiency of two. Isolated complexes 
have been excluded from the reaction graph. 

5 Conclusions 

In this paper, we have presented a computational method for determining 
weakly reversible networks with minimal deficiency which are linearly con- 
jugate to a given chemical reaction network or kinetic system. 

It was shown that, for this purpose, it is sufficient to maximize the num- 
ber of linkage classes where the linkage classes are defined to include isolated 
complexes as well as traditional linkage classes with multiple complexes. The 
proposed algorithm is based on mixed integer linear programming where the 
binary variables are used to keep track of the complexes' positions in different 
linkage classes, and the continuous variables keep track of the reaction rate 
coefficients, the conjugacy coefficients, the structure of the network, and 
the emptiness/non-emptiness of complex partitions. An additional linear 
constraint on the binary variables ensures the uniqueness of the complex 
partitioning, which dramatically improves the computational efficiency of 
the algorithm. We then applied the algorithm to several examples and were 



able to use the Deficiency Zero Theorem (Theorem 2.1) and Deficiency One 



Theorem (Theorem 2.4) to determine properties of the equilibrium set of 
the original reaction network. 

Future work in this area includes: 

1. In order to apply the algorithm outlined in this paper, it is necessary 
that the rate constants of the original reaction network be specified. 
Consequently, we may be overlooking networks with a lower deficiency 
which are linearly conjugate to a given network for certain rate con- 
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stants values but not for others. Example 1 provides a good example 
of a mechanism which permits different optimal deficiencies for differ- 
ent rate constant choices. If we were, for instance, only given the rate 



constant values specified by (22) or (23), we would not realize that 
the mechanism permitted a deficiency one conjugacy for different rate 
constant values. (Research on these 'structurally-fixed' networks was 
initiated in [22] but can thus far only solve for dynamically equivalent 
networks.) 

2. Many deficiency results do not require the networks in consideration 
to be weak reversibility. For instance, the Deficiency One Theorem 



(Theorem 2.4) and the algorithm presented in [9] only require the 
networks to have a single terminal strongly linked component within 
each linkage class. It would be useful, therefore, to adapt the algorithm 
to networks which are not necessarily weakly reversible but satisfy 
some alternative structural requirements. 
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